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SUMMARY 


Based  on  the  half-space  screen  propagator  developed  in  our  previous  project,  we  successfully 
extended  the  method  to  the  case  of  irregular  surface  topography  by  the  conformal  and 
non-conformal  topographic  transforms.  Its  validity  and  potential  applications  have  been 
numerically  demonstrated  by  comparisonsonson  with  the  boundary  element  method.  The  new 
method  can  handle  combined  effects  of  small-scale  heterogeneities  (random  media)  and  rough 
random  topography  on  Lg  wave  propagation.  It  is  also  2-3  orders  of  magnitude  faster  than  the 
finite  difference  method.  In  the  P-SV  wave  case,  reflected  plane  waves  from  a  free  surface  are 
incorporated  into  the  elastic  screen  method  for  Lg  wave  simulation.  Due  to  the  presence  of  a 
surface  wave  (Rayleigh  wave),  both  the  real  and  imaginary  parts  of  the  wave-number  must  be 
included.  Body  waves  including  reflected  and  converted  waves  are  calculated  using  real 
wavenumber  integration;  surface  waves  are  calculate  with  imaginary  wavenumber  integration. 
The  comparison  between  results  from  numerical  tests  and  the  wavenumber  integration  method 
shows  excellent  agreement.  Numerical  results  show  that  this  is  a  promising  method  for 
simulating  path  effects  in  different  regions  for  various  monitoring  discriminants  such  as  Pn/Lg  or 
Sn/Lg,  etc. 


1 


1.  INTRODUCTION 


Modeling  high-frequency  regional  wave  propagation  in  complex  crustal  waveguides  is  one  of 
the  most  challenging  problems  in  theoretical  and  computational  seismology.  A  good 
understanding  of  the  propagation,  scattering,  attenuation  and  wave-type  conversion  of  regional 
phases  is  essential  for  the  application  of  regional  waves  in  monitoring  a  CTBT  (Comprehensive 
Nuclear-Test-Ban  Treaty).  Complex  waveguide  structures,  including  rough  topography,  uneven 
Moho  discontinuities  and  volume  heterogeneities  of  all  scales,  will  affect  the  propagation  of  the 
Lg  wave.  High-frequency  regional  waves  up  to  20  Hz  have  been  observed  over  distances  ranging 
from  a  few  hundred  km  to  more  than  one  thousand  km  (e.g.,  Ni  et  ah,  1996;  Herrmann  et  ah, 
1997;  Lay  et  ah,  1999).  For  these  high-frequency  waves,  new  mechanisms  may  be  involved  in  • 
their  propagation,  scattering  and  attenuation.  The  availability  of  numerical  tools  that  can 
simulate  and  analyze  Lg  wave  propagation  under  these  environments  is  crucial.  Discrimination 
and  yield  estimation  at  regional  distances  for  very  low-yield  events  are  even  more  demanding  of 
these  tools.  To  simulate  regional  phases,  traditional  techniques  usually  use  overly  simplified  crust 
models  or  do  calculations  for  very  limited  frequencies  and  distances.  The  current  DoD/DSWA 
project  is  targeted  to  develop  a  very  flexible,  while  still  efficient,  numerical  method  for  simulating 
high-frequency,  long-distance  Lg  wave  propagation  in  a  realistic  crustal  waveguide  environment. 

[2]  Recently,  the  generalized  screen  method  has  been  introduced  into  seismic  wave 
simulations  and  applied  to  problems  of  both  exploration  seismology  and  earthquake  seismology. 
The  generalized  screen  propagator  (GSP)  is  based  on  the  one-way  wave  equation  and  the 
one-return  approximation.  The  propagator  neglects  backscattered  waves  but  correctly  handles 
all  forward  propagated  energy,  e.g.,  multiple  forward  scattering,  focusing/defocusing,  diffraction, 
interference  and  conversion  between  different  wave  types.  Significant  progress  has  been  made  in 
the  development  of  an  elastic  complex  screen  (ECS)  method  for  modeling  elastic  wave 
propagation  in  complicated  structures  (Wu,  1994,  1996;  Xie  and  Wu,  1995;  Wild  and  Hudson, 
1998;  Wu  and  Wu,  1998,  1999).  The  method  is  two  to  three  orders  of  magnitude  faster  than  the 
elastic  finite-difference  method  for  medium  size  3D  problems.  The  screen  method  has  been 
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successfully  used  for  forward  modeling  (Wu,  1994;  Wu  and  Huang,  1995;  Xie  and  Wu  1995,  1996. 
1999,  2000;  Wu  and  Wu,1999)  and  as  backpropagators  for  seismic  imaging/migration  in  both 
acoustic  and  elastic  media  (Wu  and  Xie  1994;  Huang  and  Wu  1996;  Huang  et  al.,  1999a, b;  Jin 
and  Wu,  1999;  Jin  et  al.,  1999,  Xie  and  Wu,  1998). 

In  the  crustal  waveguide  environment,  Lg  is  a  guided  wave  composed  of  upgoing  and 
downgoing  waves  bouncing  between  the  free  surface  and  other  major  geological  discontinuity'' 
such  as  the  Moho  and  Conrad.  Beyond  the  critical  angle,  they  are  consistently  dominated  by 
small  angle  waves  trapped  in  the  crustal  waveguide.  Major  wave  energy  is  carried  by  :  ward 
propagating  waves.  Therefore  the  neglect  of  backscattered  waves  in  the  propagation  will  not 
change  the  main  features  of  regional  waves  in  most  cases.  By  neglecting  backscattering  in  the 
theory,  the  method  becomes  a  forward  marching  algorithm,  the  next  step  of  propagation  only 
depends  on  the  present  values  of  the  wavefield  in  a  transverse  cross-section  and  the 
heterogeneities  between  the  two  cross-sections.  The  savings  of  computing  time  and  storage 
memory  are  enormous.  This  makes  it  a  viable  method  for  high-frequency,  long-distance  Lg  wave 
propagation  in  3D  elastic  problems. 

A  half-space  generalized  screen  propagator  has  been  introduced  (Wu  et  al.,  1996,  2000a)  to 
accommodate  the  free-surface  boundary  condition  and  to  handle  SH  wave  propagation  in 
complex  crustal  waveguides.  The  new  propagator  has  been  calibrated  extensively  against  several 
full-wave  numerical  methods  for  different  crustal  models.  In  these,  the  wavenumber  integration 
method  is  chosen  for  flat  structures  and  the  finite-difference  method  is  used  for  laterally 
heterogeneous  crustal  waveguides.  Excellent  agreement  is  obtained  and  shows  the  validity  and 
accuracy  of  the  new  one-way  method.  For  a  crust  model  with  a  propagation  distance  of  250  km 
and  a  dominant  frequency  of  0.5  Hz,  the  GSP  method  is  about  300  times  faster  than  the 
finite-difference  method  with  similar  accuracy.  The  method  also  been  used  to  simulate  Lg 
propagation  as  well  as  related  energy  partitioning  and  attenuation  in  random  media  (Wu  et  al., 
2000b).  It  is  found  that  the  leakage  attenuation  caused  by  large- angle  forward  scattering  may 
play  an  important  role  for  Lg  attenuation  and  blockage  in  some  regions.  The  equivalent  Q  of 
leakage  attenuation  versus  normalized  scale  length  ( ka )  of  random  heterogeneities  agrees  well 
with  observations  and  scattering  theory. 

In  the  current  DoD/DSWA  project,  the  half-space  GSP  method  for  SH  waves  has  been 
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extended  to  include  irregular  topography  using  both  conformal  and  non-conformal  coordinate 
transforms.  In  the  conformal  transform  method,  the  coordinate  system  is  rotated  according  to 
the  local  topographic  slope,  and  the  mirror  image  method  is  applied  to  the  local  plane  surface. 
The  non-conformal  method  is  a  surface  flattening  transform  that  transforms  the  surface 
perturbation  into  modified  volume  perturbations.  The  former  method  is  suitable  for  dealing  with 
smoothly  varying  topography,  while  the  latter  can  treat  rough  but  moderate  topography.  To 
show  their  validity  and  potential  applications,  numerical  examples  are  conducted  and  compared 
with  the  results  from  the  boundary  element  method.  A  preliminary  investigation  of  the 
combined  effects  of  small-scale  heterogeneities  (random  media)  and  rough  random  surface  on  Lg 
wave  propagation  is  carried  out. 

Another  new  development  under  this  project  is  a  2D  P-SV  elastic  complex  screen 
propagator  for  crustal  waveguides  with  a  flat  free  surface.  Plane  wave  reflection  coefficients  are 
introduced  to  calculate  down-going  waves.  Due  to  the  existance  of  a  surface  (Rayleigh)  wave, 
both  the  real  and  imaginary  parts  of  the  wavenumber  must  be  included.  Body  waves,  including 
the  reflected  and  converted  waves,  can  be  calculated  using  real  wavenumber  integration,  while 
surface  waves  can  be  calculated  with  imaginary  wavenumber  integration.  The  comparison 
between  the  results  of  the  numerical  tests  and  the  wavenumber  integration  method  shows 
excellent  agreement.  The  seismic  responses  for  a  more  realistic  crustal  model,  the  Flora- Asnes 
model,  are  simulated  using  the  elastic  screen  method. 
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2.  SH-WAVE  CASE 


A  half-space  GSP  method  for  modeling  Lg  wave  propagation  in  smoothly  varying, 
heterogeneous  crustal  waveguides  has  been  developed  by  Wu  et  al.  (2000a).  The  detailed 
derivation  of  the  method  and  numerical  tests  as  well  as  its  application  can  be  found  in  the 
papers  of  Wu  et  al.  (2000a, b).  In  this  section,  we  focus  on  extending  the  half-space  GSP  method 
to  handle  an  irregular  free  surface  through  application  of  conformal  and  non-conformal 
coordinate  transforms.  For  a  complete  description  of  the  extended  GSP  method,  the  half-space 
screen  propagator  is  briefly  summerized. 


2.1  A  BRIEF  SUMMARY  OF  THE  HALF-SPACE  SCREEN 
PROPAGATOR 

For  a  2D  SH  problem,  only  the  y-component  of  the  displacement  field,  denoted  as  u,  exists. 

With  the  perturbation  method,  the  medium  and  the  wave  field  are  decomposed  into 

p  =  Po  +  Sp\  p  =  no  +  u  =  u°  +  U 

where  po  and  p0  are  the  density  and  shear  rigidity  of  the  background  medium,  6p  and  8p  are  the 
corresponding  perturbations,  u°  is  the  primary  field  and  U  is  the  scattered  field.  The  SH  wave 
equation  in  the  frequency  domain  can  be  written  as 

p0S72u  +  u2p0u  =  —  [u>28pu  +  V  •  SpVu]  (1) 

For  each  step  of  the  marching  algorithm  under  the  forward-scattering  approximation,  the 
total  field  at  X\  is  calculated  as  the  sum  of  the  primary  field  propagating  in  the  half-space  from 
x'  to  xi,  and  the  scattered  field  caused  by  the  heterogeneities  in  the  thin-slab  between  x'  and  X\. 
The  thickness  of  the  slab  should  be  made  thin  enough  to  ensure  the  validity  of  the  local  Born 
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approximation.  The  Green’s  function  in  the  homogeneous  half-space  can  be  obtained  by  the 
image  method.  The  stress  should  vanish  at  the  free  surface  2  =  0.  We  consider  only  the  one-way 
wave  propagation  in  the  positive  x-direction.  The  scattered  field  by  the  thin-slab  can  be 
calculated  by  (Wu  et  ah,  2000a) 


U{xuKz)  =  Up(xuKz)  +  {/m(xi,A'2) 

Up(x1,I\z)  =  ik  j '  dxen(Xl~x)C[^-£p(z)u0(z)\ 

Up(xi,Kz)  =  ik dxe^x^  {C[£p{z)dxu0(z)} 

-  iS[^£fl(z)dzu0{z)]^  (2) 


where 


Po  Po 


and  7  =  \Jk?  —  is  the  propagating  wave  number  in  the  x-direction  and  I\z  is  transverse 
wavenumber  along  the  z-axis,  and 


s_iJL  s_Ii. 

‘  ikdx’  ik  dz 

In  the  above  equations,  C[f(z)\  and  <S[/(2)]  are  the  cosine  and  sine  transforms.  u0,  dxu0  and 
dzuo  can  be  calculated  by 

uo(x,z)  =  C-1[e,V^-"')u0(x,,/r')]  (3) 

8xu0(x,z)  =  C-'ie^-^uoix'Jt')} 

dzu0(x,z)  =  i<S-i[enW)^.Uo(x',A'')]  (4) 

where  C~X[F{KZ) ]  and  <S_1[F(A^)]  are  the  inverse  transforms  of  C[/(2)j  and  <S[/(2)].  The  above 
equations  are  the  general  wide-angle  formulation.  When  the  energy  of  crustal  guided  waves  are 
mainly  carried  by  small-angle  waves  (with  respect  to  the  horizontal  direction),  the  phase-screen 
approximation  can  be  invoked  to  simplify  the  theory  and  calculations.  Summing  up  the  primary 
and  scattered  fields  and  invoking  the  Rytov  transform  result  in  the  dual-domain  (phase  space), 
the  expression  of  the  phase-screen  propagator  for  this  case  is 

u(*i,  Kz)  «  e*(*i -*)c  \eikS*{z)uQ{x',  z )]  (5) 
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where  etfcSj(2)  is  the  phase  delay  operator  with 


S.(z) 


1  fxi 

2  Jx,  dxi£p(x' z)  ~  2)]  «  Axes(z). 


(6) 


Here  es(z)  is  the  average  S-wave  slowness  perturbation  over  the  thin  slab  at  depth  2, 


•w = -i—  r 

Xi  —  X1  Jx' 


with  s(x,  2)  l/n(x,<r),  and  Ax  —  (xi  —  x')  is  the  thin-slab  thickness.  Equation  (5)  is  the  SH 

phase-screen  propagator  for  the  half-space.  It  has  a  similar  form  as  the  whole  space  propagator 
with  the  Fourier  transform  replaced  by  a  cosine  transform. 


2.2  CONFORMAL  COORDINATE  TRANSFORM  FOR 
IRREGULAR  TOPOGRAPHY 

Equation  (5)  is  the  SH  phase-screen  propagator  for  the  half-space  with  a  flat  surface.  In  the  case 
of  irregular  topography,  the  global  mirror  symmetry  for  the  problem  no  longer  exists.  However, 
if  a  first  order  approximation  (local  plane-surface  approximation)  for  the  topography  is  taken,  we 
can  modify  the  mirror  image  method  to  a  local  mirror  image  method  and  apply  the 
corresponding  coordinate  transform  to  obtain  a  GSP  solution  for  this  case. 

Figure  1  shows  the  geometry  of  the  derivation.  Assume  u£  is  the  incident  field  on  S+,  then 
uo  on  S~  is  also  known  as  the  mirror  image  of  u£  about  the  local  plane  surface.  The  total 
wavefield  composed  of  S+  and  S~  is  the  sum  of  the  primary  field  which  propagates  in  the 
homogeneous  background  medium  and  the  scattered  field  which  is  generated  by  the  local 
heterogeneities  in  the  thin-slab.  The  effects  of  the  heterogeneities  and  the  topography  can  be 
calculated  separately  for  each  step  of  the  GSP  method.  The  effect  of  the  slant  free-surface  can 
be  incorporated  into  the  propagation  integral.  Assume  ut(x,z)  is  the  total  field  including  the 
scattering  effect  of  the  volume  heterogeneities.  The  field  u(xi,zi)  in  front  of  the  integral  surface 
can  be  calculated  by  the  Kirchhoff  integral 

=  js  d>  {<>(*,  *;  -  -g(-x’ „,(x, ,) j 

=  +  (7) 
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Figure  1.  Geometry  of  the  coordinate  transform. 
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where  g{.)  is  the  Green’s  function  for  the  full  space  with  homogeneous  velocity  distribution.  S  is 
the  integration  surface  and  5+  and  S~  are  the  lower  and  upper  half  surfaces,  respectively.  The 
Rayleigh  integral  can  be  used  to  replace  the  Kirchhoff  integral  for  each  half  surface  integral.  For 
the  lower  half-space  the  contribution  of  S+  is 

dg{x,z-,xi,zi) 


roo  ( 

ii+(zi,zi)  =  -2  J  dzuf(x,z)- 

=  j- JdKTeiKTZ'ut(xuKT ) 


(8) 


where 

roc 

ufixi^Kj)  =  e1'r(xi_x)  /  d2lut{x,z)e-iRTZl  (9) 

Jo 

Here  uf(x,z)  is  the  total  field  equal  to  the  sum  of  incident  field  Uq  (x,z)  and  the  scattered  field 
U+(x,z)  caused  by  the  heterogeneities  within  the  slab  between  x  and  x\  (Wu,  1994;  Wu  et  ah, 
2000a).  If  we  put  the  slab  entrance  at  x  =  x'  and  the  field  on  the  screen  S+  at  the  entrance  as 
uf(x',  z') ,  then 


ut{x',z') 

U+(x',z') 


u$(x',z')  +  U+(x',z') 

k2  [  dxe~n^l~x,)  f  dz{g(xi,zi\x,z)ep(x,z)u0(x,z) 
Jx'  Jo 

^Vg(si,  zx\  x,  z)  •  eM(®,  z)Vti0(i,  z)  j 


(10) 

(11) 


For  the  bent  upper  half  screen,  we  perform  a  coordinate  transform  by  clockwise  rotation  of 
26  to  a  new  coordinate  system  ( x,z ).  The  relation  connecting  the  two  systems  is 

x  =  x  cos  26  -F  z  sin  29  (12) 

z  =  —x  sin  26  +  z  cos  26  (13) 

In  the  new  system,  the  surface  S~  is  parallel  to  the  5-axis,  so  that 

u;(xuKt)  =  f°  d~z'u;(?,~z')e-,RT;'  (14) 

J  —  oo 

where  ui(x',-z')  =  ttt+(x',z').  The  field  in  the  space  domain  can  be  obtained  by  synthesizing 
the  contributions  from  all  the  plane  waves 

ut-(*T,z)  =  /  dKTe^  "*'>  etAV\t-  (x~',  KT)  (15) 

where 

ur(x',KT)=  1°  dz'u;(i',z’)e-ih^'.  (16) 

J  — oo 
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Transforming  back  to  the  original  coordinate  system, 

jiK  T  exp  | z  [(7  cos  26  -  I\T  sin  20) (a:  1  -  x')  +  (7  sin  26 

+I\'t  cos  20)27]  I  Ut(x':  I\'t)  (17) 

If  we  transform  the  (k?,  7)  system  into  (At,  7), 

u^(xi,zi)=  f  <IKtu'[(xi,Ktcos26  —  /ysin26)et'),(Xl~x^ethTZl  (18) 

The  total  field  is  a  summation  of  the  contributions  from  both  u+(xi,zi)  and  Ut(x\,Zi) 

u{xu2l)  =  J  dKT^(x'-x,)eiK™  [«+(*',  Kt) 

+u^(x',  Kj  cos  26  —  7  sin  20)j  (19) 

The  wavenumber  integral  can  also  be  done  by  a  FFT. 

When  small-angle  waves  prevail  such  as  in  the  case  of  Lg  propagation,  the  spectral 
interpolation  in  equation  (19),  which  is  tricky  and  even  unstable  in  some  cases,  can  be  avoided 
and  replaced  by  operations  in  the  space  domain  using  a  narrow-angle  wave  approximation.  From 
(19),  it  can  be  seen  that  to  calculate  the  reflection  response  we  need  to  find  the  spectral 
components  uf(— Kt).  We  will  try  to  obtain  the  approximate  space-domain  operations 
corresponding  to  the  wavenumber-domain  interpolation.  We  know 

POO 

u+(-Kt)=  dzei(-KTCOs2e+^sin26)zut{z)  (20) 

Jo 

With  the  narrow-angle  approximation,  7  «  k,  therefore, 

ut(—KT)  =  r  dz'eiKTZ'  [-i— e,fc{tan2fl)2'u+(27  cos  20)1  (21) 

Jo  L cos  20 

We  see  that  the  corresponding  operations  for  wavenumber-domain  interpolation  in  the 
space-domain  are  a  modulation  plus  a  coordinate  stretching. 

2.3  NON-CONFORMAL  COORDINATE  TRANSFORM  FOR 
IRREGULAR  TOPOGRAPHY 

An  alternative  to  using  the  GSP  method  for  solving  a  range- depen  dent  boundary  condition 
problem  is  to  employ  a  simple  flattening  transformation  given  by  Beillis  &  Tappert  (1979).  The 
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transformation  is  defined  as 


X  =  x 


{  C  =  2-  h(x) 

where  h(x)  is  the  height  function  of  free  surface.  Equation  (22)  shows  that  the  transformation 
gives  only  a  shift  to  depth  variable  2,  i.e.,  depth  measurement  starts  from  the  free  surface.  Under 
the  above  transformation,  the  medium  parameters  in  the  new  coordinate  {\',C}  can  be  expressed 
by  £(X,C)  =  (j-(x,z  -  h(x))  =  p(x,z)  and  /$(x,C)  =  p(z,z  ~  Kx))  ~  p{x,z)- 


Strictly  speaking,  there  is  no  one-to-one  relationship  to  the  displacement  fields  of  the  two 
different  coordinate  systems  {x,z}  and  {x,(}.  However,  if  h(x)  is  smooth  enough  and  the  height 
fluctuations  are  small,  the  displacement  field  u(x,z )  may  be  approximated  by  u(x,(),  i.e., 

u(x,z)&u(x,0  =  u(x,z-h(x))  (23) 


Substituting  eqs  (22)  and  (23)  into  the  SH  wave  equation  results  in 

C)atC)  +  cUf0  =  ~%x'  0u2ii(x,  0 

+Mx, 0  ~ 


The  simple  surface  flattening  transformation  introduces  three  extra  terms  proportional  to 
the  slope  of  the  free  surface.  At  first  glance,  it  seems  that  the  transformation  makes  the  problem 
more  complicated.  But  it  turns  a  complex  waveguide  with  an  irregular  surface  into  a  simple  one 
with  a  flat  surface.  Under  the  new  coordinate  system  {x,  C}>  the  boundary  condition  becomes 

A^(x.C  =  0)  =  0  (25) 

Using  p(x,  C)  =  Po  +  Sp(Xi  C)  and  fi(x,  C)  =  Po  +  ^(x>  C)  Equation  (24)  can  be  rewritten  as 

(v2  +  kl)  U(x,  0  =  -kl  {Mx,  0  +  Fs(x,  0}  0  (26) 

where  hi  —  u/voi  vo  —  r/po/po  and  the  operators  are 


Mx,C)  -  £p(xi()  +  pV • 

h{x’°  =  -wMh"{x)-k+2h'(x) 


cP"  cP 

wr{h'(x))2w 
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Equation  (26)  describes  SH  wave  propagation  in  a  homogeneous  medium  (background  medium), 
but  with  non-zero  equivalent  volume  forces.  The  equivalent  volume  forces  Fv(Xi  C)  and  .Fs(\sC) 
are  from  the  volume  heterogeneities  and  the  surface  flattening  transformation,  respectively.  We 
assume  that  the  change  of  the  fields  caused  by  such  volume  forces  is  small  when  compared  with 
the  primary  background  field.  Thus  perturbation  theory  can  be  applied  to  solve  equation  (26). 
Using  Green’s  theorem,  the  scattered  field  by  the  equivalent  forces  can  be  written  as 


U(Xi,h)  =  kl  f  d2Vgh(xi,kr,-,x,C) 


tV(x,0  +  Fs(x,c)]«(x,C) 


(29) 


where  the  2D  volume  integration  is  over  the  volume  V  including  all  the  heterogeneities  in  the 
modeling  space.  gh  is  a  half-space  scalar  Green’s  function  (Wu  et  ah,  2000a). 


Comparing  equations  (26)  and  (1),  we  see  that  the  only  difference  between  the  two 
equations  is  the  additional  term  Es(x,C)-  Therefore  we  can  apply  the  half-space  GSP 
(generalized  screen  propagator)  developed  by  Wu  et  al.  (2000a)  to  the  calculation  of  equation 
(26).  Under  the  forward-scattering  approximation,  for  each  step  of  the  marching  algorithm,  the 
total  field  at  Xi  is  calculated  as  the  sum  of  the  primary  field  which  is  the  field  free-propagated  in 
the  half-space  from  x'  to  Xu  and  the  scattered  field  caused  by  the  heterogeneities  in  the 

J 

thin-slab  between  x!  and  Xi-  The  thickness  of  the  slab  should  be  made  thin  enough  to  ensure  the 
validity  of  the  local  Born  approximation.  Under  this  condition,  the  Green’s  function  can  be 
approximated  by  the  homogeneous  half-space  Green’s  function.  Using  the  image  method,  Wu  et 
al.  (2000a)  derived  a  half-space  Green’s  function  as 


9oiXi,k6X,0  =  ^-en(xi  x)2cos(kc() 
27 


(30) 


The  scattered  field  at  the  exit  xi  of  the  thin  slab  by  the  volume  heterogeneities  Fv  is 


Uty  (X..*t)  =  '-~£‘  ***<*'-*'  {C 


-1,(0 MO  ~  4(C)8*6o(C) 

7 


V 


+&  ^(O^o(C)j  |  (31) 

where  uq(Q  is  the  incident  field  at  the  entrance  x'  of  the  thin  slab.  In  a  similar  way,  we  can 
obtain  the  scattered  field  at  the  exit  Xi  by  the  equivalent  force  Fs 

Uh(x„k)  =  Y  /)'  dxertn-tic  {  ^  [-2h\x)A  +  (h\x)?B  -  h"(x)D ] )  (32) 
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where 


A  =  a~wr =  ~iS~'  (33) 

b  =  — =  -c-1  fcvv'>->'iao(x-it0]  (34) 

D  =  =  _5-i  [fe;e'VI»-«')Ao(x^')j  (35) 

Equations  (31-35)  are  the  dual-domain  expressions  of  the  wide-angle  thin-slab  propagator  for  SH 
problems  in  an  arbitrary  half-space  model  with  arbitrary  heterogeneities  and  an  irregular  surface. 

Small-angle  approximation 

When  the  energy  of  crustal  guided  waves  is  carried  mainly  by  small-angle  waves  (with 
respect  to  the  horizontal  direction),  the  small-angle  approximations  can  be  invoked  to  simplify 
the  theory  and  calculations.  For  small-angle  waves,  ^  <  7  w  )'  k  ko ,  Upv(xAc)  can  be 
approximated  by 

VK  (x,  h)  «  ikoehAxC  lS.(OMx\  0]  (36) 

where  Ay  =  Xi  ~  X1  is  thickness  of  the  thin  slab,  and 

s,(0  =  \  dx  fc(x,  0  -  M  x,  0)  (37) 

For  Ups(x-i^)i  substituting  equations  (33-  35)  into  (32),  we  have 

U?s(xuk()  =  ir  f“dXe^-x)  r  r  dk[e‘^-^eik'<( 

.  *  Z‘ y  Jo  fJ 0  J— 00 

j^2h'(x)fc^7/  —  ih"(x)h ^  —  h\x)2k (-2j  tio(x,>  k<)  (38) 

For  a  sufficiently  smooth  h{x)  with  respect  to  x,  under  the  small-angle  approximation,  the  last 
two  terms  in  the  above  equation  are  high-order  small  quantities.  Neglecting  the  last  two  terms, 
we  obtain 

UpAxuh)  =  l~enAx  f  d(2cos(k<Q^-  [  dx2h'(x)  f  dk[ekAk[u0{x ,  k'c) 

s  l  Jo  fj  0  Jx'  J-OO 

=  -2(x1)e"Axc|^5-1[*:;4„(x',lcJ)]|  (39) 

with 

2(xi)  =  M*i)  ~  h{x')  (40) 
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It  is  clear  that  under  the  small-angle  approximation  the  scattered  field  Up  is  propotional  to  the 
height  difference  of  the  adjacent  two  screens  for  each  forward  step.  For  an  upward  slope, 

Z(x)  <  0.  the  scattered  field  is  in-phase  with  the  background  field  and  strengthens  the 
background  field,  while  for  a  downward  slope,  Z(x)  >  0,  the  scattered  field  is  out-phase  with  the 
background  field  and  weakens  the  background  field.  The  total  field  can  be  obtained  through 
summing  the  primary  field  which  freely  propagates  in  the  background  medium  and  the  scattered 
fields  (36)  and  (39) 

u(Xuk)  =  e,'|,AxC  |fi0(y,,C)  +  ik0Ss(()u0(x',()  —  [fcju0(x\*c)]} 

l  /*  o  J 

«  e^xcje'^IOC-1  [ua(x',y()]-  Z(xi)^-S-'  [*'«„(*',*')]  j  .  (41) 

Equation  (41)  is  the  dual-domain  expression  of  the  small-angle  approximation  for  SH  wave 
simulation  in  complex  waveguides  with  arbitrary  heterogeneities  and  irregular  surfaces. 

2.4  NUMERICAL  VERIFICATIONS  AND  SIMULATION 
EXAMPLES 

The  validity  and  accuracy  of  the  GSP  method  applied  to  uniform  crustal  waveguides  have  been 
tested  by  comparing  to  the  finite  difference  and  reflectivity  methods  (Wu,  et.  ah,  2000a, b).  In 
this  section  we  show  some  examples  demonstrating  the  validity  and  the  potential  of  the  GSP 
with  coordinate  transforms  applied  to  crustal  waveguides  with  topography  compared  to  the 
boundary  element  method. 

Figure  2  shows  the  synthetic  seismograms  for  a  Gaussian  hill  model  using  the  conformal 

(g—gp)2 

screen  method.  The  Gaussian  hill  is  represented  by  h(x)  =  —h0e  2 J  with  x0  =  62.25 km, 
h0  =  4 km,  a  =  9.129 km.  Synthetic  seismograms  calculated  with  the  more  accurate  boundary 
element  method  (Fu  and  Wu,  1999)  are  also  given  as  a  reference  .  The  solid  lines  represent  the 
screen  method  and  the  dashed  lines  represent  the  boundary  element  method.  The  comparison 
indicates  that  the  screen  method  gives  a  satisfactory  result.  It  correctly  modeled  waveforms 
between  60  and  70  km  distance,  where  two  reflections  from  the  convex  free  surface  interfere  with 
each  other  and  generate  complex  waveforms.  Note  that  the  coordinate  stretch  z/cos26  increases 
very  fast  with  large  angle  6 ,  the  conformal  screen  method  works  only  for  smoothly  varying 
topography. 
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Horizontal  distance  (km) 


Figure  2.  Synthetic  seismograms  using  the  conformal  transform  method  (solid  lines).  For  the  calculation, 
dx=dz=0.25km  and  dt=0.05sec.,  the  source  is  located  at  a  horizontal  location  of  0  km,  a  depth 
of  32km  and  the  dominant  frequency  of  source  time  function  is  3Hz.  Receivers  are  on  the  free 
surface.  The  dashed  lines  are  calculated  with  boundary  element  (BE)  method  and  used  as  a 
reference  (They  overlie  the  solid  lines  very  closely).  For  the  calculation  of  the  BE  method, 
each  wavelength  contains  5  boundary  elements  at  least  (Fu  and  Wu,  1999). 
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Figure  3  shows  the  synthetic  seismograms  using  the  non-conformal  screen  method  for  the 
Gaussian  hill  model  used  in  Figure  2.  The  excellent  agreement  between  the  screen  method  and 
boundary  element  method  is  clear  except  in  the  vicinity  of  the  hill  top  where  a  small  discrepancy 
exists  both  in  wave  shapes  and  amplitudes.  The  error  will  decrease  when  the  step  length  Ax 
reduces.  For  forward  marching  algorithms,  the  step  length  Ax  may  be  adjusted  according  to  the 
roughness  of  topography.  The  more  severe  the  topography,  the  finer  the  Ax  should  be.  The 
non-conformal  screen  method  can  handle  more  severe  topography  than  the  conformal  screen 
method. 

Figure  5  shows  a  comparison  of  synthetic  seismograms  and  attenuation  between  the  screen 
method  and  the  boundary  element  method  applied  to  a  crustal  waveguide  with  a  rough  random 
surface  shown  in  Figure  4.  For  the  rough  random  surface,  the  correlation  length  is  2.5km,  RMS 
height  fluctuation  is  0.6km  (range  -1.5  to  1.75  km),  the  mean  height  is  -0.11  km.  Figure  5  shows 
a  comparison  of  synthetic  seismograms  calculated  by  the  non-conformal  screen  method  and  the 
BE  method,  and  the  corresponding  energy  attenuation  curves,  respectively.  The  source  is  located 
at  a  depth  of  8  km,  the  dominant  frequency  of  the  source  time  function  is  1Hz.  The  thick 
smoothly  varying  curve  in  Figure  5b  is  the  energy  distribution  for  a  flat  free  surface.  We  see  that 
except  for  large- angle  Moho  reflections,  the  results  of  the  screen  method  agree  well  with  those  of 
the  BE  method.  However,  for  this  example,  the  screen  method  took  about  35  minutes,  but  the 
BE  method  took  about  72  hours. 

With  the  extension  of  the  screen  propagator  for  handling  surface  topography,  Lg  wave 
propagation  in  more  realistic  crustal  waveguides  with  both  volume  heterogeneities  and  surface 
topography  can  be  simulated.  Figures  6-8  show  an  investigation  of  the  combined  effects  of  rough 
topography  and  volume  heterogeneity  on  Lg  wave  propagation  using  the  extended  screen 
method.  The  rough  topography  is  the  same  as  shown  in  Figure  4.  The  heterogeneities  are  only 
for  velocity  variation.  The  heterogentity  correlation  lengths  are  6  km  in  range  and  4  km  in  depth, 
RMS  are  5%  and  10%,  respectively.  The  source  is  located  at  a  depth  of  8  km.  The  dominant 
frequencies  of  the  source  time  function  are  0.5Hz  in  Figure  6,  1Hz  in  Figure  7  and  2Hz  in 
Figure  8,  respectively.  The  corresponding  dominant  wavelengths  are  7  km,  3.5  km  and  1.75  km, 
respectively.  The  thickly  dashed  line  calculated  by  the  finite-difference  method  for  an  uniform 
crustal  waveguide,  is  used  as  a  reference.  We  see  from  Figures  6-8  that  random  heterogeneities 
combined  with  rough  topography  drastically  increase  the  attenuation  of  high  frequency  Lg  waves. 
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Horizontal  distance  (km) 


Figure  3.  Synthetic  seismograms  using  the  non-conformal  transform  method  (solid  lines).  For  the 
calculation,  dx=0. 125km,  dz=0.25km  and  dt=0.05sec.  Source  is  located  at  the  depth  of  32km 
and  the  dominant  frequency  of  source  time  function  is  3Hz.  Receivers  are  on  free  surface.  The 
dashed  lines  are  calculated  with  the  boundary  element  (BE)  method  and  used  as  a  reference. 
For  the  calculation  of  the  BE  method,  each  wavelength  contains  5  boundary  elements  at  least 
(Fu  and  Wu,  1999). 
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Figure  4.  A  crustal  model  with  rough  random  surface.  The  correlation  length  is  2.5km,  RMS  height 
fluctuation  is  0.6  km  (range  -1.5  to  1.75  km),  the  mean  height  is  -0.11  km. 


Attenuation  and  variation  of  Lg  waves  with  distance  are  related  to  frequency,  the  characteristic 
length  of  random  heterogeneities,  and  roughness  of  the  surface.  This  example  shows  that  the 
screen  method  can  handle  the  effects  of  both  volume  heterogeneities  and  moderately  rough 
topography  on  Lg  wave  propagation  at  long  propagation  distances  and  high  frequencies. 


2.5  BOUNDARY  ELEMENT  METHOD  FOR  REAL  CRUSTAL 
WAVEGUIDES 

During  this  project,  we  also  studied  Lg  wave  attenuation  due  to  rough  surfaces  in  two 
representative  geological  profiles  in  Western  China  using  the  boundary  element  (BE)  method. 
The  topographic  and  Moho  depth  data  for  these  waveguides  are  from  Fan  and  Lay  (1998).  The 
first  path  is  from  the  epicenter  of  June  29,  1995  event  (focal  depth  15  km,  latitude  51.923N  and 
longitude  103.075E)  to  station  WMQ.  The  second  path  is  from  the  epicenter  of  June  2,  1990 
event  (focal  depth  17  km,  latitude  32.45N  and  longitude  92.81E)  to  station  WMQ.  These 
waveguide  models  are  shown  in  Figures  9  and  10,  respectively.  The  source  time  function  used  for 
calculating  synthetic  seismograms  is  a  Ricker  wavelet  with  a  dominant  frequency  of  1.0  Hz.  For 
the  second  event,  Lg  blockage  has  been  observed,  while  for  the  first  event,  no  Lg  blockage  is 
observed.  From  Figure  9,  we  see  that  the  Moho  discontinuity  of  the  first  waveguide  is  relatively 
flat.  The  crustal  thickness  is  larger  (over  40  km)  and  rather  uniform  throughout  the  entire 
waveguide.  The  surface  roughness  is  locally  dramatic  (over  2  km)  but  generally  moderate.  The 
energy  decay  curve  in  Figure  9  shows  strong  fluctuations  correlated  with  the  topography.  The 
overall  shape  of  the  energy  decay  curve  is  relatively  stable,  especially  beyond  1000  km  where  the 
energy  has  been  trapped  in  the  waveguide  and  the  attenuation  is  weak.  The  model  with  Lg 
blockage  shows  a  completely  different  waveguide  structure.  The  large  lateral  variations  of 
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Figure  5.  Comparison  between  the  non-conformal  screen  method  and  BE  method  for  a  crustal  waveg¬ 
uide  with  a  rough  random  surface,  (a)  Synthetic  seismograms,  (b)  Energy  distribution  with 
horizontal  distance. 
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Figure  6.  Lg  wave  attenuation  versus  horizontal  distances.  A  random  medium  whose  correlation  lengths 
are  6km  in  range  and  4 km  in  depth,  and  RMS  velocity  fluctuations  are  5%  and  10%,  re¬ 
spectively.  Source  is  located  at  a  depth  of  8km,  the  dominant  frequency  (fO)  of  source  time 
function  is  0.5Hz.  ’’rough”  in  the  Figure  means  the  crust  with  rough  topography,  and  ”ho” 
and  ”het”  mean  the  crusts  which  are  homogeneous  and  inhomogeneous,  respectively. 


Figure  7.  Lg  wave  attenuation  versus  horizontal  distances.  A  random  medium  whose  correlation  lengths 
are  6 km  in  range  and  4 km  in  depth,  and  RMS  velocity  fluctuations  are  5%  and  10%,  re¬ 
spectively.  Source  is  located  at  a  depth  of  8 km,  the  dominant  frequency  (fO)  of  source  time 
function  is  1Hz.  ’’rough”  in  the  Figure  means  the  crust  with  rough  topography,  and  ”ho”  and 
”het”  mean  the  crusts  which  are  homogeneous  and  inhomogeneous,  respectively. 
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Figure  8.  Lg  wave  attenuation  versus  horizontal  distances.  A  random  medium  whose  correlation  lengths 
are  6 km  in  range  and  4 km  in  depth,  and  RMS  velocity  fluctuations  are  5%  and  10%,  re¬ 
spectively.  Source  is  located  at  a  depth  of  8 km,  the  dominant  frequency  (fO)  of  source  time 
function  is  2Hz.  ’’rough”  in  the  Figure  means  the  crust  with  rough  topography,  and  ”ho”  and 
”het”  mean  the  crusts  which  are  homogeneous  and  inhomogeneous,  respectively. 

topography  and  the  shallowing  Moho  discontinuity  rapidly  change  the  crust  into  a  necking 
structure  after  800  km  (Figure  10).  The  anomalous  crustal  necking  may  be  responsible  for  the 
Lg  blockage.  Although  current  calculations  are  limited  to  lower  frequencies  (less  than  2  Hz),  the 
energy  decay  curve  shows  a  rapid  attenuation  near  the  necking  structure. 
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A  sample  profile  from  epicenter  to  WMQ  of  the  June  29,  1995  enevt. 


Energy  attenuation  versus  distance 


Figure  9.  Energy  attenuation  versus  distance  for  the  path  without  Lg  wave  blockage  in  the  Tibet 
region.  The  model  is  from  the  event  of  June  29,  1995,  with  the  location  of  51.923N 
and  103.075E,  the  magnitude  of  5.6  and  the  source  depth  of  15  km. 
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A  sample  profile  from  epicenterto  WMQ  of  the  June  2, 1990  event. 
Distanee(km) 


Figure  10.  Energy  attenuation  versus  distance  for  the  path  without  Lg  wave  blockage  in  the  Tibet 
region.  The  model  is  from  the  event  of  June  2,  1990,  with  the  location  of  32.45N  and 
92.81E,  the  magnitude  of  5.5  and  the  source  depth  of  17  km. 
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3.  P-SV  CASE 


Another  new  development  in  this  project  is  that  the  2D  P-SV  elastic  complex  screen 
method  has  been  applied  to  Lg  wave  simulation  in  crustal  waveguides  with  a  flat  surface  by 
incorporating  plane  wave  reflection  coefficients  in  the  elastic  complex  screen  method.  Body 
waves  including  the  reflected  and  converted  waves  can  be  calculated  by  real  wavenumber 
integration,  surface  waves  (Rayleigh  waves)  can  be  done  with  imaginary  wavenumber  integration. 
Numerical  tests  show  excellent  agreement  with  the  theory. 

For  P-SV  elastic  wave  one-way  propagation  in  heterogeneous  crustal  waveguides,  the 
derivation  and  application  of  screen  propagators  is  much  more  complicated.  Unlike  the  case  of 
SH  waves,  the  image  method  of  generating  the  half-space  Green’s  function  cannot  be  used  to 
account  for  the  effect  of  free  surface. 

Figure  11  illustrates  the  principle  of  the  method.  Since  the  screen  method  is  a  forward 
marching  algorithm,  the  re-application  of  reflection  coefficients  for  each  forward  step  will 
contribute  extra  energy,  already  taken  in  the  last  step,  to  reflections.  To  avoid  extra  energy 
contributions  to  reflections,  we  extend  the  model  twice  in  the  vertical  direction  from  the  free 
surface;  the  extended  part  has  the  parameters  of  the  background  medium  and  will  be  used  to 
keep  records  of  the  upgoing  waves  which  can  be  used  for  the  calculation  of  reflected/converted 
waves  by  the  free  surface.  The  complex  screen  method  will  be  applied  to  such  an  extended 
model  for  elastic  wave  simulation.  For  each  forward  step  we  apply  the  reflection  coefficient 
formulas  to  the  incident  field  (only  to  upgoing  waves)  to  calculate  the  reflections  and  conversions 
from  the  free  surface.  Figure  1  also  shows  such  a  process  by  an  upgoing  ray  from  source. 
Knowing  the  reflections  and  conversions  ,  we  calculate  the  scattered  field  of  the  reflected  fields 
using  the  complex  screen  method  and  then  combine  the  scattered  field  into  the  incident  field  to 
obtain  the  new  incident  field.  In  this  section,  a  simple  description  of  the  complex  screen  method 
and  reflection  coefficient  formulas  is  given. 
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Figure  11.  Illustration  of  the  screen  method. 


3.1  COMPLEX  SCREEN  PROPAGATOR  FOR  ELASTIC  WAVES 


A  complete  derivation  of  the  complex  screen  method  for  modeling  elastic  waves  can  be  found  in 
Wu  (1994).  A  short  review  of  the  method  for  forward  scattering  is  given  in  this  section.  Suppose 
that  the  parameters  of  the  elastic  medium  and  the  total  wave  field  can  be  expressed  as 


p(x)  =  p0  +  Sp(y.) 
a(x)  =  a0  +  5a(x) 
/3(x)  =  fa  +  50{x) 

u(x)  =  U0(x)  +  U(x) 


where  p0,  e*o  and  0O  are  the  parameters  of  the  background  medium,  <5p(x),  <5a(x)  and  80{x)  are 
the  corresponding  perturbations,  u0(x)  and  U(x)  are  the  incident  field  and  the  scattered  field, 

X  =  xex  +  zez  is  a  2-D  position  vector  in  Cartesian  coordinates.  We  take  ex  to  be  the  primary 
propagation  direction  perpendicular  to  the  thin  slab.  The  incident  field  Uo(x)  at  x\  e.g.,  the 
entrance  of  the  thin  slab,  can  be  decomposed  into  a  superposition  of  plane  P  and  S  waves 

u0(x)  =  ~ ^  J  dk  [u£(fc,a:')  +  Uo(fc,z')]  etkz  (42) 

where  k  is  the  transverse  wavenumber  of  plane  waves.  Superscripts  P  and  S  denote  P  and  S 
waves,  respectively.  The  forward  propagated  field  is  composed  of  the  primary  wave  and  forward 
scattered  P  and  S  waves.  Therefore,  at  Xi  the  exit  of  the  slab,  it  can  be  expressed  as 


\  f 

u/(z,a:i)  =  —  J  dk  [u^fc,  a^)  +  usf(k,  xx) 


Jkz 


(43) 
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where 

u -/(Mi)  =  {e,'7oAru^(fc,a;,)  +  u^p(A;,x1)  +  ufp(fc,a:1)}  (44) 

u  /(Mi)  =  {ei^Axu^{k,x')  +  usfs(k,x1)  +  ups(k,x1)}  (45) 

where  and  ^p  are  x  components  of  P  and  S  wavenumbers  in  the  background  medium. 

Ax  =  x\  —  x'  is  the  thickness  of  the  thin  slab.  Subscript  /  denotes  forward  scattering,  and  PP , 
PS,  SP  and  SS  indicate  the  scattering  between  different  wave  types.  Under  the  complex  screen 
approximation,  equations  (44)  and  (45)  can  be  expressed  by 

uPj{k,x  0  =  e^A*up(M0  +  uf(Mi) 

=  kehaAx  I  dzeikzup{x',  z)e-^(x',z)AX  +  (46) 


where 


where 


and 


u j (k,  x^) 


ups{k,x  x) 

uSsP{k,Xl) 


e^Axus0(k,x')  +  uPS(k,x  1) 

en0Axkp  x  i^kp  x  J  dze~tkz  x  (2,  x/)e_’fc^(a7'’2)Aa:| 

+ups(k,Xl) 

ikpe'lf>Ax  Axr]{Ax)kp  x  x  J  dzelkz  x  up(z,x') 

Po  _  l\  &p(z,x')  +  2Po  SP(z,x')'  | 

Qo  2  y  po  c*o  /?o  J 

-i^e’^Ax^Ax)*  •  j  dze~ikzui(z,x') 

^0  _  l\  fy(z,x')  +  2(3q5(3{z,x') 

«o  2  /  po  do  ft 0 


rj(Ax)  =  sine 


{kp  fca) 


Ax 


-  k{kp-ka)Ax 


e  2 


(47) 


(48) 


(49) 


(50) 


7.  =  (*J  -  fc2)‘/2,  7«s  =  (*i  -  P)‘/2  •  (51) 

=  w/qo  and  kp  =  u;//30  are  P  and  S  wavenumbers,  a0  and  /30  are  P  and  S  velocities  of  the 
background  medium.  The  unit  wavenumber  vectors  ka  and  kp  are  given  by 


koc  —  k ot/k(x ,  kp  —  \&.pjkp 


(52) 
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Equations  (46-49)  are  the  dual-domain  formulation  of  the  complex  screen  method.  The  incident 
fields  Uq  and  Uq  interact  with  screens  (i.e.,  heterogeneities)  in  the  space  domain,  but  propagate 
in  the  wavenumber  domain.  We  combine  the  common  type  scatterings  (P-P  and  5-5)  with 
incident  fields  (/up  and  Uq)  into  pure  phase  screens  to  ensure  stability  of  the  method  for  models 
with  large  contrasts. 

3.2  FREE-SURFACE  REFLECTIONS  OF  PROPAGATING  P-SV 
WAVES 

From  the  screen  method,  the  incident  P  and  5  waves  at  x  =  x'  can  be  expressed  as  a 
superposition  of  plane  waves  by 

/OO 

dkuZ(k,x')eikz  (53) 

-OO 

/OO 

dkn^(k,x')e'kz.  (54) 

-OO 

Applying  the  reflection  coefficients,  the  reflected  P  and  S  waves  due  to  an  incident  P  wave  can 
be  expressed  by 

u pp(x,z)  =  eh°{x-x>)  J  dA>£(AT,^l^aie~'A'iZ  (55) 

ups(x,z)  =  e*ba(*-*')  J  dKz\u%{Kz,x')\PS*2e-iK*z  (56) 

where  la  =  \fk2  —  K\  is  the  propagating  wavenumber  for  P  waves  (here  in  the  x-direction)  and 
k;  =  jk}  -ki  +  ki  is  the  transverse  wavenumber  of  the  converted  S  wave  determined  by 
Snell’s  law.  Unit  vectors  are  a4  =  (7 Q,—Kz)/kQ  and  a2  =  ( K*,ia)/kp .  |up (Kz,x')\  is  the  scalar 
spectrum  of  incident  P  wave  with  transverse  wavenumber  Kz  (here  in  the  z-direction).  The 
reflected  P  and  S  waves  due  to  the  incident  plane  S  wave  can  be  obtained  by 

uSP(x,z)  =  e^(*-*')  j dKz\ul{Kz,x')\SPkze-iK''z  (57) 

uss(x,z)  =  J  dAT|u£(AT,x')j55a4e-tp*2  (58) 

where  ip  =  \Jk20~  I\2Z  is  the  propagating  wavenumber  for  S  waves  (here  in  the  x-direction)  and 
K'*  =  yjk2  —  kp  +  K2  is  the  transverse  wavenumber  of  the  reflected  P  wave.  Unit  vectors  are 
»3  =  ilfh -K'*)/ka  and  a4  =  ( Kz ,  l/3)/kp.  |Uq(AT,x')|  is  the  scalar  spectrum  of  incident  S  wave 

/  N  /  \  /  X  /V 

with  transverse  wavenumber  Kz.  PP ,  PS,  SP  and  55  in  equations  (55—58)  are  reflection 
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Figure  12.  Reflection  coefficients  (in  log)  on  free  surface  versus  horizontal  slowness.  Elastic  half-space 
has  a  =  5km/ s  and  f3  =  3km/ s.  pa  corresponds  to  P  slowness  and  ps  to  S  slowness. 


coefficients  for  the  free  surface  (Aki  and  Richards,  1980).  Figure  12  is  an  example  of  those 
reflection  coefficients  versus  horizontal  slowness  (ray  parameter  p).  In  Figure  12  pa  corresponds 
to  P  slowness  (inverse  velocity)  and  ps  to  S  slowness.  For  p  <  pa,  P  and  S  waves  are  both 
homogeneous  waves,  their  transverse  wavenumbers  are  real.  For  p  >  pb,  P  and  S  waves  are  both 
inhomogeneous  waves,  their  transverse  wavenumbers  are  imaginary.  For  Pa  <  P  <  Pb,  P  wave  is 
inhomogeneous  while  S  wave  is  homogeneous.  A  Rayleigh  pole  is  located  at  the  region  of  p  >  pb- 
In  general,  we  can  calculate  all  of  the  reflected  waves  using  equations  (55—58),  once  the  incident 
fields  |uq  |  and  |u„  |  are  known.  However,  numerically,  it  is  more  convenient  to  separate  the 
calculation  of  equations  (55—58)  into  homogeneous  and  inhomogeneous  waves,  respectively. 

For  homogeneous  waves,  equations  (55)  and  (58)  (common-type)  can  be  implemented  by  the 
FFT.  However,  the  reflected  waves  of  converted-type  cannot  be  directly  implemented  by  the 
FFT  because  a  nonlinear  relationship  exists  between  Kz  and  K*  for  P-S  (or  Kz  and  K'*  for  S-P). 
Although  we  can  obtain  uniform  samples  with  respect  to  Kz  and  K*  (or  Kz  and  K'*)  by 
complex  variable  interpolation  for  application  of  the  FFT,  numerical  tests  have  shown  that  the 
noise  due  to  the  interpolation  is  so  strong  that  the  accumulated  errors  increase  very  fast  with 
multiple  step  propagation.  In  our  study,  the  direct  summations  over  the  incident  waves  (p  <  pa 
for  P  incidence  or  p  <  Pb  for  S  incidence)  are  performed  for  calculating  the  converted  reflections. 
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3.3  SURFACE  WAVE  PROPAGATION  USING  SCREEN 


PROPAGATOR 

For  inhomogeneous  waves  such  as  the  Rayleigh  wave,  transverse  wavenumbers  are  imaginary  so 
that  equations  (55—58)  cannot  be  calculated  by  FFT.  However,  the  imaginary  transverse 
wavenumber  makes  the  propagation  of  inhomogeneous  waves  simple.  The  phase  advance  is  only 
along  the  horizontal  direction.  It  is  easily  incorporated  into  the  screen  method  in  the 
wavenumber  domain,  once  the  spectra  of  the  inhomogeneous  waves  are  known.  Another 
important  feature  of  the  propagation  of  inhomogeneous  waves  is  the  exponential  decay  occuring 
only  in  the  direction  perpendicular  to  propagation.  Then  the  spectra  of  inhomogeneous  waves 
can  be  calculated  with  the  Laplace  transform. 


3.4  NUMERICAL  TESTS  AND  EXAMPLES 

Figure  13  shows  synthetic  seismograms  calculated  with  equations  (55—58)  for  an  elastic 
half-space  using  only  homogeneous  waves.  A  point  explosion  source  is  located  at  a  depth  of  16 
km  and  the  dominant  frequency  of  source  time  function  is  1Hz.  The  first  4  receivers  are  placed 
along  the  free  surface  separated  from  the  source  by  100~124  km,  and  the  last  5  receivers  are 
placed  in  vertical  profile  at  132  km  with  depths  ranging  from  0  to  32  km.  The  results  calculated 
with  the  wavenumber  integration  (WI)  method  (dashed)  are  also  shown  as  references.  Since  the 
source  is  deep  compared  with  the  propagation  distance,  the  Rayleigh  wave  is  very  weak  in  the 
exact  solution.  Figure  13a  shows  the  vertical  component  of  the  displacement  and  13b  shows  the 
horizontal  component.  From  Figure  13  we  see  that  the  calculations  of  the  reflection  and 
conversion  by  the  free  surface  are  in  excellent  agreement  with  the  theory. 

Figure  14  shows  snapshots  generated  by  the  screen  method  for  the  Flora-Asnes  crustal 
model  in  the  NORSAR  region.  The  parameters  for  this  model  are  listed  in  Table  1.  The  model 
has  a  low  velocity  layer  in  the  top  1  km  and  a  velocity  discontinuity  at  depth  15  km.  The 
receivers  are  on  the  surface.  A  double-couple  source  is  located  at  the  depth  of  16km  and  has  a 
dominant  frequency  of  2Hz  (Ricker  wavelet).  We  see  that  both  P  and  S  waves  are  well  excited. 
Most  events  (direct  P  and  S,  head  waves,  PmP,  PmS,  pS,  mantle  wave  etc.)  can  be  clearly 
identified.  Figure  15  is  the  corresponding  seismograms. 
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Vertical  displacement 


Time  (s) 


Figure  13.  Synthetic  seismograms  calculated  by  the  elastic  screen  method  (solid)  and  wavenumber 
integration  method  (dashed)  for  an  elastic  halfspace.  Only  homogeneous  waves  are  involved 
in  the  results  of  elastic  screen  method,  (a)  shows  the  vertical  components  of  displacement,  (b) 
shows  the  horizontal  components.  A  point  explosion  source  is  located  at  the  depth  of  16km 
and  the  dominant  frequency  of  source  time  function  is  1Hz.  The  first  4  receivers  are  placed 
along  the  free  surface  separated  from  the  source  by  100~124km,  and  the  last  5  receivers  are 
placed  in  vertical  profile  at  132km  with  depths  ranging  0~32km. 
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Table  1.  Flora-Asnes  crust  model 


Thickness  (km) 

a(km/s) 

0(km/s) 

p(g/cm3) 

1.00 

5.20 

3.46 

2.80 

6.51 

3.76 

Infinity 

8.05 

4.65 

3.30 

Figure  16  shows  a  similar  result  to  Figure  14  for  the  Flora-Asnes  crustal  model,  but  an 
explosive  source  is  used  to  instead  of  the  double-couple  source  and  the  source  is  located  at  a 
depth  of  2  km.  Figure  17  is  the  corresponding  seismograms.  We  see  that  the  seismic  responses 
for  these  two  different  sources  are  fairly  different,  especially  for  the  Lg  parts. 

Figures  18  shows  an  example  of  such  a  treatment  for  Rayleigh  wave  propagating  in  a 
homogeneous  elastic  half-space.  The  source  is  located  at  a  depth  of  2  km  and  has  a  dominant 
frequency  of  0.5  Hz.  The  vertical  receiver  array  is  at  a  horizontal  distance  of  100  km.  (a)  shows 
the  vertical  component  and  (b)  shows  the  horizontal  component  of  Rayleigh  waves.  The  dotted 
lines  are  the  exact  solution.  The  agreement  between  the  screen  calculation  and  the  theory  in 
excellent. 


31 


■200  km- 


I-*-  -  ' 


Snapshots  for  Flora-Asnes  crustal 
model.  A  double  couple  source  is 
located  at  a  depth  of  16  km 


Figure  14.  Snapshots  (horizontal  component  of  displacement)  for  Flora-Asnes  crustal  model  using 
P-SV  elastic  screen  method.  A  double-couple  source  is  located  at  a  depth  of  16  km  and 
has  a  dominant  frequency  of  2  Hz.  The  thicknesses  of  layers  (from  top  to  bottom)  are  1 
km,  14  km,  22  km  and  infinity,  respectively. 
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Synthetic  seismograms  for  Flora-Asnes  crust  model 

(Double-couple  source,  16km  depth) 


Figure  15.  Synthetic  seismograms  for  Flora-Asnes  crust  model  (Table  1).  A  double-couple  source 
is  located  at  a  depth  of  16 km  and  the  source  time  function  has  a  dominant  frequency 
of  1Hz 
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Snapshots  for  Flora-Asnes  crust  model 

(An  explosive  source,  2km  depth) 


200km 


Figure  16.  Snapshots  for  Flora-Asnes  crust  model  (Table  1).  An  explosive  source  is  located  at  a 
depth  of  2 km  and  the  source  time  function  has  a  dominant  frequency  of  2 Hz 
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Figure  17.  Synthetic  seismograms  for  Flora-Asnes  crust  model  (Table  1). 
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Figure  18.  Comparison  of  synthesized  Rayleigh  wave  (dotted)  using  screen  method  with  the  exact 
solution  (solid).  Source  is  located  at  the  depth  of  2km  and  has  a  dominant  frequency  of 
0.5Hz.  (a)  shows  the  horizontal  components  of  displacement  of  Rayleigh  wave  and  (b)  shows 
the  vertical  components.  The  half-space  parameters  are  a  =  6km/ s  and  (3  =  3.5 km/s. 
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4.  CONCLUSION 


Based  on  the  half-space  screen  propagator  developed  in  our  previous  project,  we  successfully 
extended  the  method  to  the  case  of  irregular  surface  topography  by  conformal  or  non-conformal 
topographic  transforms.  The  validity  and  potential  of  the  extended  method  has  been  numerically 
demonstrated  by  comparing  the  results  with  the  boundary  element  method.  The  method  can 
handle  the  combined  effects  of  small-scale  heterogeneities  (random  media)  and  rough  random 
topography  on  Lg  wave  propagation.  For  medium  size  problems,  the  screen-propagator  method 
is  2-3  orders  of  magnitude  faster  than  the  finite  difference  method. 

In  the  case  of  P-SV  elastic  screen  propagators,  plane  wave  reflection  calculations  have  been 
incorporated  into  the  elastic  screen  method  for  Lg  wave  simulation.  Due  to  the  presence  of  the 
surface  wave  (Rayleigh  wave),  all  wavenumber  components  (real  or  imaginary  parts)  must  be 
included.  Body  waves  including  the  reflected  and  converted  waves  can  be  calculated  by  real 
wavenumber  integration,  surface  waves  (Rayleigh  waves)  can  be  done  with  imaginary 
wavenumber  integration.  Numerical  tests  against  the  wavenumber  integration  method  show 
excellent  agreement.  An  example  of  simulated  Lg  wave  propagation  in  Flora- Asnes  crust  model 
has  demonstrated  the  promise  of  simulating  path  effects  in  different  regions  for  various 
discriminants  in  the  monitoring  system,  such  as  P„/Ls,  Sn/Lg,  etc. 
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ATTN:  R.  MORROW,  ROOM  5741 

US  GEOLOGICAL  SURVEY 
905  NATION  CENTER 
12201  SUNRISE  VALLEY  DR 
RESTON,  VA  20192 

ATTN:  W.  LEITH 

US  GEOLOGICAL  SURVEY 
345  MIDDLEFIELD  RD  MS  977 
MENLO  PARK,  CA  94025 

ATTN:  S.  DETWEILER 
ATTN:  DR.  W.  MOONEY 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

BATTELLE 

MANAGER,  ENERGETIC  SYSTEMS  & 
SECURITY  TECHNOLOGIES 
505  KING  AVE 

COLUMBUS,  OH  43201-2693 

ATTN:  NEAL  OWENS  (7-2-081) 

BBN  CORPORATION 

1300  N  17TH  STREET,  SUITE  400 

ARLINGTON,  VA  22209 

ATTN:  DR.  D.  NORRIS  ' 

ATTN:  R.  GIBSON 
ATTN:  J.  PULLI 

CENTER  FOR  MONITORING  RESEARCH 
1953  GALLOWS  ROAD,  SUITE  260 
VIENNA,  VA  221 82  3997 

ATTN:  DR.  K.  L.  MCLAUGLHLIN 
ATTN:  DR.  R.  WOODWARD 
ATTN:  DR.  R.  NORTH 
ATTN:  DR.  X.  YANG 
ATTN:  LIBRARIAN 

ENSCO,  INC. 

5400  PORT  ROYAL  ROAD 
P.O.  BOX  1346 

SPRINGFIELD,  VA  22151  2312 
ATTN:  D.  BAUMGARDT 
ATTN:  Z.  DER 

WESTON  GEOPHYSICAL  CORPORATION 
27  BEDFORD  ST,  SUITE  102 
LEXINGTON,  MA  02420 

ATTN:  DR.  D.  REITER 
ATTN:  J.  LEWKOWICZ 
ATTN:  DR.  A.  ROSCA 
ATTN:  DR.  I.  TIBULEAC 
ATTN:  M.  JOHNSON 


ITT  INDUSTRIES 
ITT  SYSTEMS  CORPORATION 
1680  TEXAS  STREET  SE 
KIRTLAND  AFB,  NM  87117  5669 
2  CYS  ATTN:  DTRIAC 
ATTN:  DARE 

TITAN  CORPORATION  (ATS) 

1900  CAMPUS  COMMONS  DR  SUITE  600 
RESTON,  VA  20191-1535 

ATTN:  DR.  C.  P.  KNOWLES 

MISSION  RESEARCH  CORPORATION 
8560  CINDERBED  ROAD,  SUITE  700 
NEWINGTON,  VA  22122 

ATTN:  DR.  M.  FISK 

MULTIMAX,  INC 

1441  MC  CORMICK  DRIVE 

LANDOVER,  MD  20785 

ATTN:  DR.  I.  N.  GUPTA 
ATTN:  W.  RIVERS 

MULTIMAX,  INC 

1090  N  HIGHWAY  A1A  SUITE  D 

INDIANLATIC,  FL  32903 

ATTN:  DR.'H.  GHALIB 

SCIENCE  APPLICATIONS  INTERNATIONAL 
CORP 

10260  CAMPUS  POINT  DRIVE 
SAN  DIEGO,  CA  92121  1578 
ATTN:  DR.M.  ENEVA 
ATTN:  DR.  G.  E.  BAKER 
ATTN:  DR.  J.  STEVENS 
ATTN:  DR.  D.  ADAMS 

SCIENCE  APPLICATIONS  INT’L  CORP 
1227  S.  PATRICK  DR  SUITE  110 
SATELLITE  BEACH,  FL  32937 
ATTN:  DR.  M.  FELIX 
ATTN:  DR.  H.  GIVEN 


URS  CORPORATION 
566  EL  DORADO  STREET 
PASADENA,  CA  91109  3245 
ATTN:  DR.  C.  SAIKIA 
ATTN:  DR.  G.  ICHINOSE 

DIRECTORY  OF  OTHER  (LIBRARIES  AND 
UNIVERSITIES) 

BOSTON  COLLEGE 
INSTITUTE  FOR  SPACE  RESEARCH 
140  COMMONWEALTH  AVENUE 
CHESTNUT  HILL,  MA  02167 

ATTN:  DR.  D.  HARKRIDER 
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BROWN  UNIVERSITY 

DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
75  WATERMAN  STREET 
PROVIDENCE,  Rl  02912  1846 

ATTN:  PROF.  D.  FORSYTH 

CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
DIVISION  OF  GEOLOGY  &  PLANETARY 
SCIENCES 

PASADENA,  CA  91 125 

ATTN:  PROF.  DONALD  V. 
HELMBERGER 

ATTN:  PROF.  THOMAS  AHRENS 

UNIVERSITY  OF  CALIFORNIA  BERKELEY 
281  MCCONE  HALL 
BERKELY,  CA  94720  2599 

ATTN:  PROF.  B.  ROMANOWICZ 
ATTN:  DR.  D.  DREGER 

UNIVERSITY  OF  CALIFORNIA 
DEPT  OF  STATISTICS 
DAVIS,  CA  95616 

ATTN:  R.H.  SHUMWAY,  DIV 
STATISTICS 

UNIVERSITY  OF  CALIFORNIA  SAN  DIEGO 
SCRIPPS  INSTITUTION  OF  TECHNOLOGY 
9500  GILMAN  DRIVE 
LA  JOLLA,  CA  92093  0225 

ATTN:  DR.L.  DEGROOT  -  HEDLIN 
ATTN:  DR.  M.  HEDLIN 
ATTN:  PROF.  F.  VERNON 
ATTN:  PROF.  J.  BERGER 

UNIVERSITY  OF  CALIFORNIA  SANTA  CRUZ 
INSTITUTE  OF  TECTONICS 
SANTA  CRUZ,  CA  95064 

ATTN:  DR.  R.  S.  WU 
ATTN:  PROF.  T.  LAY 

UNIVERSITY  OF  COLORADO  BOULDER 
DEPT  OF  PHYSICS,  CAMPUS  BOX  390 
BOULDER,  CO  80309 

ATTN:  DR.R.  ENGDAHL 
ATTN:  M.  RITZWOLLER 
ATTN:  PROF.  C.  ARCHAMBEAU 

COLUMBIA  UNIVERSITY 
LAMONT  DOHERTY  EARTH  OBSERVATORY 
PALISADES,  NY  10964 
ATTN:  DR.  J.  XIE 
ATTN:  DR.  W.  Y.  KIM 
ATTN:  PROF.  P.G.  RICHARDS 
ATTN:  DR.  M.  TOLSTOY 

UNIVERSITY  OF  CONNECTICUT 
DEPARTMENT  OF  GEOLOGY  &  GEOPHYSICS 
STOORS,  CT  06269  2045 

ATTN:  PROF.  V.F.  CORMIER,  U-45, 
ROOM  207 


CORNELL  UNIVERSITY 

DEPARTMENT  OF  GEOLOGICAL  SCIENCES 

3126  SNEE  HALL 

ITHACA,  NY  14853 

ATTN:  PROF.  M.A.  BARAZANGI 

HARVARD  UNIVERSITY 
HOFFMAN  LABORATORY 
20  OXFORD  STREET 
CAMBRIDGE,  MA  02138 

ATTN:  PROF.  A.  DZIEWONSKI 
ATTN:  PROF.  G.  EKSTROM 

INDIANA  UNIVERSITY 

DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
1005  10TH  STREET 
BLOOMINGTON,  IN  47405 

ATTN:  PROF.  G.  PAVLIS 

IRIS 

1200  NEW  YORK  AVENUE,  NW  SUITE  800 
WASHINGTON,  DC  20005 

ATTN:  DR.  D.  SIMPSON 

IRIS 

1408  NE45TH  ST  #201 
SEATTLE,  WA  98105 

ATTN:  DR.  T.  AHERN 

MASSACHUSETTS  INSTITUTE  OF 
TECHNOLOGY 

EARTH  RESOURCES  LABORATORY 
42  CARLETON  STREET 
CAMBRIDGE,  MA  02142 

ATTN:  DR.  W.  RODI 
ATTN:  PROF.  M.N.  TOKSOZ 

MICHIGAN  STATE  UNIVERSITY  LIBRARY 
450  ADMINISTRATION  BUILDING 
EAST  LANSING,  Ml  48824 
ATTN:  K.  FUJITA 

NEW  MEXICO  STATE  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
LAS  CRUCES,  NM  88003 

ATTN:  PROF.  J.  Nl 
ATTN:  PROF.  T.  HEARN 

NORTHWESTERN  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
1847  SHERIDAN  RD 
EVANSTON,  IL  60208 

ATTN:  PROF.  E.  OKAL 

PENNSYLVANIA  STATE  UNIVERSITY 
GEOSCIENCES  DEPARTMENT 
403  DEIKE  BUILDING 
UNIVERSITY  PARK,  PA  16802 

ATTN:  PROF.  C.  AMMON 
ATTN:  PROF.  S.  ALEXANDER 
ATTN:  DR.  A.  NYBLADE 
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SAN  DIEGO  STATE  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
SAN  DIEGO,  CA  92182 

ATTN:  PROF.  S.  M.  DAY 


SOUTHERN  METHODIST  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
P.O.  BOX  750395 
DALLAS,  TX  75275 

ATTN:  B.  STUMP 
ATTN:  E.  HERRIN 
ATTN:  P.  GOLDEN 

UNIVERSITY  OF  HAWAII-  MANOA 
P.O.  BOX  1599 
KAILUA-KONA.HI  96745  1599 

ATTN:  DR.  M.  A.  GARCES 

UNIVERSITY  OF  MISSISSIPPI 
1  COLISEUM  DRIVE 
UNIVERSITY,  MS  38677 

ATTN:  PROF.  H.  BASS 

UNIVERSITY  OF  SOUTHERN  CALIFORNIA 
520  SEAVER  SCIENCE  CENTER 
UNIVERSITY  PARK 
LOS  ANGELES,  CA  90089  0483 

ATTN:  PROF.  C.  G.  SAMMIS 
ATTN:  PROF.  T.  JORDAN 

UNIVERSITY  OF  WISCONSIN  MADISON 
1215W  DAYTON  ST 
MADISON,  Wl  53706  1600 

ATTN:  DR.  C.  THURBER 

ST  LOUIS  UNIVERSITY 
EARTH  &  ATMOSPHERIC  SCIENCES 
STATION  3507  LACLEDE  AVE 
ST  LOUIS,  MO  63103 

ATTN:  PROF.  B.  J.  MITCHELL 
ATTN:  PROF.  R.  HERRMAN 

UNIVERSITY  OF  MEMPHIS 
3904  CENTRAL  AVE 
MEMPHIS,  TN  38152 

ATTN:  DR.  J.  PUJOL 

UNIVERSITY  OF  MEMPHIS 
3876  CENTRAL  AVE 
MEMPHIS,  TN  38152 

ATTN:  DR.  C.  LANGSTON 


UNIVERSITY  OF  TEXAS  EL  PASO 
DEPT  OF  GEOLOGICAL  SCIENCES 
EL  PASO,  TX  79901 

ATTN:  PROF.  G.  KELLER 
ATTN:  DR.  D.  DOSER 
ATTN:  DR.  A.  VELASCO 

FOREIGN 

AUSTRALIAN  GEOLOGICAL  SURVEY 
ORGANIZATION 

CORNER  OF  JERRAGOMRRRA  & 
NINDMARSH  DRIVE 
CANBERRA,  ACT  2609 
AUSTRALIA 

ATTN:  D.  JESPSON 

GEOPHYSICAL  INSTITUTE  OF  ISRAEL 
POB  182 

LOD,  7100  ISRAEL 

ATTN:  DR.  Y.  GITTERMAN 
ATTN:  DR.  A.  SHAPIRA 

GEOLOGICAL  SURVEY  OF  CANADA 
7  OBSERVATORY  CRESCENT 
OTTAWA  K1A  0Y3  ONT 
CANADA 

ATTN:  C.  WOODGOLD 

I.R.I.G.M.  B.P.  68 

38402  ST.  MARTIN  D’HERES 

CEDEX,  FRANCE 

ATTN:  DR.  M.  BOUCHON 

MINISTRY  OF  DEFENSE 
PROCUREMENT  EXECUTIVE 
BLACKNESS,  BRIMPTON 
READING  FG7-4RS  ENGLAND 

ATTN:  DR.  D.  BOWERS 

NTNF/NORSAR 
P.O.  BOX  51 

N-2007  KJELLER,  NORWAY 

ATTN:  DR.  F.  RINGDAL 
ATTN:  T.  KVAERNA 
ATTN:  S.  MYKKELTVEIT 


UNIVERSITY  OF  TEXAS  AUSTIN 
IGS130 

AUSTIN,  TX  78712 

ATTN:  DR.  J.  PULLIAM 

UNIVERSITY  OF  TEXAS  AUSTIN 
IGS131 

AUSTIN,  TX  78712 

ATTN:  DR.  M.  SEN 
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